
# Packages ----------------------------------------------------------------

# using a reproducible environment
renv::restore()

# pacman allows to check/install/load packages with a single call
# if (!require("pacman")) install.packages("pacman") # already in renv.lock
library("pacman")

# packages to load (and install if needed) -------------------------------
pacman::p_load(
  here,      # easy file paths
  see,       # theme_modern and okabeito palette
  report,    # reporting various info 
  labelled,  # labelled data
  # ---- packages specific to this project ----
  scales,
  lme4,
  car,
  simr,
  readxl,
  openxlsx,
  statmod,
  qqplotr,
  emmeans,
  latex2exp,
  ggbeeswarm,
  ggpubr,
  patchwork,
  quarto,
  easystats,
  # ---
  tidyverse  # modern R ecosystem
)


# Custom functions shared across scripts ----------------------------------
source(here("scripts/_functions.R"))


# Global cosmetic theme ---------------------------------------------------

theme_set(theme_modern(base_size = 14)) # from see in easystats

# setting my favourite palettes as ggplot2 defaults
options( 
  ggplot2.discrete.colour   = scale_colour_okabeito,
  ggplot2.discrete.fill     = scale_fill_okabeito,
  ggplot2.continuous.colour = scale_colour_viridis_c,
  ggplot2.continuous.fill   = scale_fill_viridis_c
)


# Fixing a seed for reproducibility ---------------------------------------
set.seed(14051998)


# Adding all packages' citations to a .bib --------------------------------
knitr::write_bib(c(.packages()), file = here("bibliography/packages.bib"))

1 Rationale

We hypothesised that difference in reaction times (RT) between congruent and incongruent trials (the congruence effect) in the tasks would be modulated by the visual imagery ability of participants, and therefore that we would observe between-group differences in the magnitude of this effect. In a statistical model with categorical predictors, evidence in favour of our hypothesis would be provided by a significant two-way interaction between Group and Congruence condition. Estimating the statistical power to capture such an effect, given an experimental design, a sample size and an effect size, requires simulations from a model close to the one that will be used for the analysis of real data, as well as assumptions about the structure of data. The power estimation procedure requires 3 steps:

Specifying the structure and coefficients of a generative model
Generating a synthetic dataset using the generative model and an expected effect size
Fitting the same model on the simulated data and test the statistical significance of the interaction of interest

Steps 2 and 3 are then repeated as many times as necessary to reach a desired precision of estimation. The proportion of synthetic datasets for which the effect of interest is significant is an estimation of power. Steps 2 and 3 can also be repeated for a range of effect sizes, if the expected effect size is not known exactly which is what we have done here. We detail and justify below every assumption made for the generative model.

2 Generative model

In real datasets, the distributions of RTs are right-skewed and can be modelled with transformations using Gamma or inverse Gaussian distributions. However, such a modelling choice is only an approximation as RTs do not follow exactly any distribution of the exponential family for which generalized models can be fitted. It is difficult to evaluate the impact of such approximation on type I and type II errors. Thus, for power calculation, we generated normally distributed RT and modelled them with a regular mixed model, which simplifies the model parametrization as well as alleviates the computational burden.

The generative model has been created using the simr R package.

Creating a dataset following the structure of our data for the generative model

# First creating a dataframe with no response data, but whose structure reflects the dataset that will be collected. 

# vectors for subjects
id <- seq(1:150)
groups <- rep(c("aphantasia", "control"), 75)

# conditions
congruences <- list("congruent", "incongruent")
colors <- list("colored", "uncolored")
# repeated 4 times each
congruence = rep(congruence, 4)
color = rep(colors, 4)

# creating the simulation df
df <- 
  tibble(
    id = id, 
    group = groups,
    congruence = list(congruence),
    color = list(color)
    ) |> 
  # crossing conditions for 64 trials/subject total
  unnest_longer(congruence) |> 
  unnest_longer(color)

Fixed effects: We included the Group (aphantasic, control), Congruence condition (congruent or incongruent) and Color condition (color or uncolored) along with all their two and three way interactions as fixed categorical predictors. Regression coefficients were set in such a way as to produce a pattern of means roughly following our hypotheses, where the control group would be faster solely in the congruent condition.

simr requires regression coefficients in order to simulate data. The functions below transform a pattern of cell means to regressions coefficients of a linear model, using a simple trick: modelling a dataframe of hypothesized cell means, as if each of them was one noise-less observation of a combination of experimental factors, provides the desired coefficient estimates.

Creating regression coefficients through models

# creating dataframes of regression coefficients
extract_coefficients <- function(effectsize) {
  df_temp <- tribble(
          ~group,   ~congruence,       ~color,   ~y,
    "aphantasia",   "congruent",    "colored",  -30,
    "aphantasia",   "congruent",  "uncolored",    0,
    "aphantasia", "incongruent",    "colored",  -30, 
    "aphantasia", "incongruent",  "uncolored",    0,
       "control",   "congruent",    "colored", -30 - effectsize,
       "control",   "congruent",  "uncolored",   0 - effectsize,
       "control", "incongruent",    "colored",  -30,
       "control", "incongruent",  "uncolored",    0
  )
  
  lm_coef <- lm(y ~ (group + congruence + color)^3, data = df_temp)

  # Set numerical errors to 0
  coef(lm_coef)[abs(coef(lm_coef)) < 10^(-10)] <- 0
  
  # converting to ms upon return
  return(coef(lm_coef)/1000)
}

Random intercepts: We included a random intercept for participants with a standard deviation of 150ms. This number was derived from the literature (see Fucci et al., 2023) - however, previous simulation studies have shown that the magnitude of random intercepts has no effect on statistical power whatsoever (Brysbaert & Stevens, n.d.), a result that we could verify in our own simulations.
Random slopes: We included random by-participant slopes for Congruence and Color, thus specifying the full model given our factors and resulting in the most conservative power estimates. Standard deviations of random slopes are virtually never reported in scientific publications; we set them at 50ms, which is in the same range as the maximum expected population-average effect size.
Correlations between random effects were all set to 0, but an alternative value (0.5) produced similar results.

This random effect structure has been specified and provided to a simr modelling function to simulate the models that will be fitted on our data.

Specifying the random effect structure and the final generative model

## Fixed effects: initialized at 0
fixed <- c(700,   # intercept
           0,0,0, # main effects
           0,0,0, # 2-way interactions
           0     # 3-way interactions
           )/1000

## Random intercepts for participants
v1 <- 0.150

## Random intercept, slopes (main effects only) and their covariances
slope  <- .05
correl <- 0
v2 <- matrix(
  c(1, correl, correl,
    correl, 1, correl,
    correl, correl, 1),
  3)

v2 <- cor_to_cov(v2, sd = c((v1), slope, slope))

## residual std dev
res <- .15

model <- makeLmer(
  formula = y ~ (group + congruence + color)^3 + (congruence + color|id),
  fixef   = fixed, 
  VarCorr = list(v2), 
  sigma = res,
  data  = df
  )

3 Simulations

Effect sizes: For the interaction Group x Congruence, the effect size can be defined as the double difference (aphantasics - controls) x (congruent - incongruent). Assuming, based on previous works, that the strongest difference would not exceed 50-60ms, and that differences below 10ms would hardly be meaningful, we ran simulations for differences between 10 and 60ms, by steps of 5ms. As the power reached the ceiling of 100% after 40ms effect sizes, we report here the results up to this size.
Number of simulations: For each effect size, we ran as many simulations as were needed to obtain a confidence interval below ±10 points around the mean. Going by steps of 50, we eventually ran 200 simulations for each effect size. These simulations ran for more than 8 hours, therefore their results are stored in .RDS R objects for ease of access: data/r-data-structures/power-analyses-list.RDS and data/r-data-structures/power-analyses-table.RDS. They do not run automatically when rendering the .qmd file.

Running the simulations for various effect sizes

# ─── Running simulations with parallel processing ─────────────────────────────

# finding the available cores for parallel processing
n_cores <- parallel::detectCores() - 1

# creating the cluster of cores
parallel_cluster <- 
  parallel::makeCluster(
    n_cores,
    type = "PSOCK"
  )

# registering the cluster for `foreach`
doParallel::registerDoParallel(cl = parallel_cluster)
# checking
# foreach::getDoParRegistered()
# foreach::getDoParWorkers()

# initializing, ready for takeoff
nsim <- 200

# Run calculations for a range of effect sizes
(simulations <-
  foreach(
    effectsize = c(10, 15, 20, 25, 30, 35, 40), 
    .packages = c("tidyverse", "simr")
    ) %dopar% {
      # Calculate regression coefficients corresponding to the effect size and 
      # update the model to simulate from accordingly
      fixef(model) <- extract_coefficients(effectsize)
      
      # Simulate data and estimate power
      powerSim(
        model, 
        nsim = nsim,
        test = simr::fixed(
          xname = "group:congruence", 
          method = 'anova'
          ))
    }) |> system.time() -> time_simulating

# stopping the cluster when we're done
parallel::stopCluster(cl = parallel_cluster)

4 Results

The analysis yielded a power of 82% starting at effect sizes around 25ms, going up to 91% at 30ms and 97% at 35ms. Every effect size equal and above 40ms reached a ceiling of 100% power. The complete results are summarized in Table S4.1 and Figure S4.1.

Generating the summary table for the power analysis

simulations_table <- 
  readRDS(file = "data/r-data-structures/power-analyses-table.RDS")

simulations_table |> 
  # All the code below is simple formatting of the table that is in the RDS
  rename(
    "Effect size (difference in ms)" = effectsize,
    "Successes"   = successes,
    "Simulations" = trials,
    "Power" = mean
  ) |> 
  mutate(across(c(lower:upper), ~round(.x, digits = 2))) |> 
  unite("CI", lower:upper, sep = ", ") |> 
  mutate(CI = paste("[", CI, "]")) |> 
  rename("95% CI" = CI) |> 
  display()

Table S4.1: Results of the power analysis through simulation.

Effect size (difference in ms)	Successes	Simulations	Power	95% CI
10	46	200	0.23	[ 0.17, 0.29 ]
15	78	200	0.39	[ 0.32, 0.46 ]
20	136	200	0.68	[ 0.61, 0.74 ]
25	164	200	0.82	[ 0.76, 0.87 ]
30	182	200	0.91	[ 0.86, 0.95 ]
35	195	200	0.97	[ 0.94, 0.99 ]
40	200	200	1.00	[ 0.98, 1 ]

Generating the plot of the power analysis results

simulations_table |> 
  ggplot(aes(
    x = effectsize,
    y = mean,
    ymin  = lower,
    ymax  = upper,
    label = effectsize
    )) +
  # the band representing the 95% CIs
  geom_ribbon(fill = "aquamarine", alpha = .2) +
  # the means and CIs at each specific effect size
  geom_pointrange(colour = "aquamarine4") +
  geom_point(colour = "aquamarine4") +
  # line connecting the dots
  geom_line(colour = "aquamarine4") +
  # horizontal line at 80% power and corresponding annotation
  geom_hline(yintercept = .8, color = "black", linetype = 1) +
  annotate(
    geom  = "text",
    label = "80% power",
    fontface = "italic",
    color = "black",
    x = -Inf,
    y = .8,
    vjust = -.5,
    hjust = -.3,
    size  = 4 
    ) +
  # horizontal dotted line at 90% power and corresponding annotation
  geom_hline(yintercept = .9, color = "gray60", linetype = 2) +
  annotate(
    geom  = "text",
    label = "90% power",
    fontface = "italic",
    color = "grey60",
    x = -Inf,
    y = .9,
    vjust = -.5,
    hjust = -.3,
    size  = 4 
    ) +
  # Aesthetics
  scale_x_continuous(
    name = "Effect size (RT difference in ms)",
    breaks = unique(simulations_table$effectsize)
  ) +
  scale_y_continuous(
    name = "Statistical power",
    breaks = seq(0, 1, .1),
    limits = c(0, 1)
    ) +
  theme_modern() +
  theme(panel.grid.major = element_line())

Figure S4.1: Results of the power analysis using simulations. 200 simulations have been computed for each effect size.

Session information for reproducibility

═════════════════════════════════════════════════════════════════════════

Analyses were conducted using the R Statistical language (version 4.4.1; R Core
Team, 2024) on Windows 11 x64 (build 22631)

Packages used:

  - quarto (version 1.4.4; Allaire J, Dervieux C, 2024)
  - qqplotr (version 0.0.6; Almeida A et al., 2018)
  - lme4 (version 1.1.35.5; Bates D et al., 2015)
  - Matrix (version 1.7.0; Bates D et al., 2024)
  - effectsize (version 0.8.9; Ben-Shachar MS et al., 2020)
  - ggbeeswarm (version 0.7.2; Clarke E et al., 2023)
  - car (version 3.1.2; Fox J, Weisberg S, 2019)
  - carData (version 3.0.5; Fox J et al., 2022)
  - statmod (version 1.5.0; Giner G, Smyth GK, 2016)
  - simr (version 1.0.7; Green P, MacLeod CJ, 2016)
  - lubridate (version 1.9.3; Grolemund G, Wickham H, 2011)
  - ggpubr (version 0.6.0; Kassambara A, 2023)
  - labelled (version 2.13.0; Larmarange J, 2024)
  - emmeans (version 1.10.3; Lenth R, 2024)
  - parameters (version 0.22.0; Lüdecke D et al., 2020)
  - performance (version 0.12.0; Lüdecke D et al., 2021)
  - easystats (version 0.7.2; Lüdecke D et al., 2022)
  - see (version 0.8.4; Lüdecke D et al., 2021)
  - insight (version 0.20.1; Lüdecke D et al., 2019)
  - bayestestR (version 0.13.2; Makowski D et al., 2019)
  - modelbased (version 0.8.8; Makowski D et al., 2020)
  - report (version 0.5.8; Makowski D et al., 2023)
  - correlation (version 0.8.5; Makowski D et al., 2022)
  - latex2exp (version 0.9.6; Meschiari S, 2022)
  - here (version 1.0.1; Müller K, 2020)
  - tibble (version 3.2.1; Müller K, Wickham H, 2023)
  - datawizard (version 0.11.0; Patil I et al., 2022)
  - patchwork (version 1.2.0; Pedersen T, 2024)
  - R (version 4.4.1; R Core Team, 2024)
  - pacman (version 0.5.1; Rinker TW, Kurkiewicz D, 2018)
  - openxlsx (version 4.2.5.2; Schauberger P, Walker A, 2023)
  - ggplot2 (version 3.5.1; Wickham H, 2016)
  - forcats (version 1.0.0; Wickham H, 2023)
  - stringr (version 1.5.1; Wickham H, 2023)
  - tidyverse (version 2.0.0; Wickham H et al., 2019)
  - readxl (version 1.4.3; Wickham H, Bryan J, 2023)
  - dplyr (version 1.1.4; Wickham H et al., 2023)
  - purrr (version 1.0.2; Wickham H, Henry L, 2023)
  - readr (version 2.1.5; Wickham H et al., 2024)
  - scales (version 1.3.0; Wickham H et al., 2023)
  - tidyr (version 1.3.1; Wickham H et al., 2024)

═════════════════════════════════════════════════════════════════════════

References

Brysbaert, M., & Stevens, M. (n.d.). Power analysis and effect size in mixed effects models: A tutorial. Journal of Cognition, 1(1), 9. https://doi.org/10.5334/joc.10

Fucci, E., Abdoun, O., Baquedano, C., & Lutz, A. (2023). Ready to help, no matter what you did: Responsibility attribution does not influence compassion in expert Buddhist practitioners. Preprint. https://doi.org/10.31234/osf.io/d6y4w